(in-package :parameterscout)

;; FIXME: should sensitivity perhaps be measured as a percentage of
;; the variable values rather than an absolute measure? Percentage
;; measure allows for the fact that resistor tolerances are measured
;; in percentages, but absolute sensitivity is needed to take into
;; account external sources of error which do not depend on the
;; variable values. Ideally, both would be offered -- but this
;; threatens to severely complicate the interface.

(defmacro sensitivity-vector-form (expr variables solution-vector)
  "Given an expression, a list of variables, and a solution
vector in the same format as those returned by the solver code,
with variable bindings in the solution vector given in the same
order as the variable list, return a vector of sensitivities. The
sensitivity vector will be one element shorter than the solution
vector, since it will not include the first, extraneous
element. The solution vector can be any form which evaluates to
the necessary value."
  (once-only (solution-vector)
    `(let ,(loop for variable in variables 
	      for n from 1
	      collect `(,variable (elt ,solution-vector ,n)))
       (declare (optimize (speed 3) (safety 0) (debug 0))
		(type double-float ,@variables))
       (vector
	,@(loop for variable in variables collect
	       (differentiator:diff expr variable))))))

(defun overall-sensitivity (vec)
  "Given a vector of n sensitivities, calculate the overall
  sensitvity. If you think of the sensitivity vector as being a
  gradient vector in n-space, then this function calculates the
  magnitude of that gradient vector. It's hoped that by choosing
  the smallest gradient vector we can get the most robustness."
  (sqrt (reduce #'+ (map-into vec #'(lambda (x) (* x x)) vec))))

(sensitivity-vector-form (/ r1 r2)
			 (r1 r2)
			 #(2.2d0 22.0d0 10.0d0))

(sensitivity-vector-form (/ (* r1 r2) (+ r1 r2))
			 (r1 r2)
			 #(1228.9156626506024d0 6800.0d0 1500.0d0))

(disassemble #'(lambda (vec)
		 (sensitivity-vector-form (/ (* r1 r2) (+ r1 r2))
					  (r1 r2)
					  vec)))